library(RTMB)
dbl_logistic_sel <- function(age, p1, p2, p3) {
gamma1 <- p1 + p2
gamma2 <- 2 * p1 + p2 + p3
asc <- 1 / (1 + exp(-log(19) * (age - gamma1) / p1))
desc <- 1 - 1 / (1 + exp(-log(19) * (age - gamma2) / p3))
(asc * desc * 0.95^(-2))
}
lognorm_from_median_cv <- function(median, cv) {
sigma <- sqrt(log(cv^2 + 1))
mu <- log(median)
list(mu = mu, sigma = sigma)
}
age <- seq(1, 15, by = 0.25)
sel_obs <- c(
0, 0, 0, 0, 0, 0, 0, 0.01, 0.01, 0.037142857, 0.064285714, 0.091428571,
0.118571429, 0.145714286, 0.172857143, 0.2, 0.205882353, 0.211764706,
0.217647059, 0.223529412, 0.229411765, 0.235294118, 0.241176471,
0.247058824, 0.252941176, 0.258823529, 0.264705882, 0.270588235,
0.276470588, 0.282352941, 0.288235294, 0.294117647, 0.3, 0.31, 0.32,
0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 0.41, 0.42, 0.43,
0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.5125, 0.525, 0.5375, 0.55,
0.5625, 0.575, 0.5875, 0.6, 0.6125, 0.625, 0.6375, 0.65, 0.6625,
0.675, 0.6875, 0.7, 0.7125, 0.725, 0.7375, 0.75, 0.7625, 0.775,
0.7875, 0.8, 0.807894737, 0.815789474, 0.823684211, 0.831578947,
0.839473684, 0.847368421, 0.855263158, 0.863157895, 0.871052632,
0.878947368, 0.886842105, 0.894736842, 0.902631579, 0.910526316,
0.918421053, 0.926315789, 0.934210526, 0.942105263, 0.95, 0.954545455,
0.959090909, 0.963636364, 0.968181818, 0.972727273, 0.977272727,
0.981818182, 0.986363636, 0.990909091, 0.995454545, 1, 1, 1, 1, 1, 1,
1, 0.989285714, 0.978571429, 0.967857143, 0.957142857, 0.946428571,
0.935714286, 0.925, 0.914285714, 0.903571429, 0.892857143, 0.882142857,
0.871428571, 0.860714286, 0.85, 0.839285714, 0.828571429, 0.817857143,
0.807142857, 0.796428571, 0.785714286, 0.775, 0.764285714, 0.753571429,
0.742857143, 0.732142857, 0.721428571, 0.710714286, 0.7
)
sel_sd <- rep(0.05, length(age))
p1_prior <- lognorm_from_median_cv(median = 1.5, cv = 0.4)
p3_prior <- lognorm_from_median_cv(median = 2.5, cv = 0.6)
p2_prior <- list(mu = 4, sigma = 1)
make_prior_data <- function(p1_prior, p2_prior, p3_prior) {
list(
age = age,
sel_obs = sel_obs,
sel_sd = sel_sd,
mu_log_p1 = p1_prior$mu,
sd_log_p1 = p1_prior$sigma,
mu_p2 = p2_prior$mu,
sd_p2 = p2_prior$sigma,
mu_log_p3 = p3_prior$mu,
sd_log_p3 = p3_prior$sigma
)
}
data <- make_prior_data(p1_prior, p2_prior, p3_prior)
parameters <- list(
log_p1 = log(1.5),
p2 = 4,
log_p3 = log(2.5)
)
make_obj <- function(data, parameters) {
f <- function(parms) {
getAll(data, parms, warn = FALSE)
p1 <- exp(log_p1)
p3 <- exp(log_p3)
sel_hat <- dbl_logistic_sel(age, p1, p2, p3)
nll <- 0
nll <- nll - dnorm(log_p1, mu_log_p1, sd_log_p1, log = TRUE)
nll <- nll - dnorm(p2, mu_p2, sd_p2, log = TRUE)
nll <- nll - dnorm(log_p3, mu_log_p3, sd_log_p3, log = TRUE)
ADREPORT(c(p1 = p1, p2 = p2, p3 = p3))
nll
}
MakeADFun(f, parameters)
}
obj <- make_obj(data, parameters)
opt <- nlminb(obj$par, obj$fn, obj$gr)